-- This example gives formulas for quadratic and cubic roots and uses the image command to visualize a complex function

------------------- quadratic roots ---------------
-- Algorithm for real parts of roots is from by W.H. Press, S. Teukolsky et al,"Numerical Recipes".
-- Occasionally 4ac << b, so one of the roots is (erroneously) called 0.
-- This formulation avoids the problem.
-- implemented by David Derbes (derbes@uhuru.uchicago.edu) for MathPad

-- Given a quadratic of the form

--          a*x^2 + b*x + c = 0
 
-- with real coefficients, find the (possibly complex) roots.
-- Roots are given in the form x + iy.
~
sgn(x) = 1 when x >= 0,-1 otherwise
D = (b*b - 4*a*c)	-- discriminant
x1 = -(b + sgn(b)*sqrt(D))/(2*a) when D >= 0, -b/(2*a) otherwise
x2 = c/x1 when D >= 0, -b/(2*a) otherwise
y1 = 0 when D >= 0, sqrt(-D)/(2*a) otherwise
y2 = -y1
~

------------------- cubic roots -------------------
-- Tartaglia's & Cardano's formulae for the roots of a cubic
-- from Universal Encyclopaedia of Mathematics
-- implemented by David Derbes for MathPad, 8 Sept 1993
-- Given a cubic of the form
--
--            a0*x^3 + a1*x^2 + a2*x + a3 = 0
-- 
-- with real coefficients, find the (possibly complex) roots.

c1 = a1/a0  -- "normalize", i.e. make leading coefficient = 1
c2 = a2/a0
c3 = a3/a0

-- discriminant D; if D > 0, one real root; else three real roots
-- (at most two distinct real if D = 0)

a = c2 - (c1*c1)/3.0
b = (((2.0*c1*c1*c1) - (9.0*c1*c2))/27.0) + c3
D = (b*b/4.0) + (a*a*a/27.0)

numRealRoots = 3 when D < 0, 1 otherwise

dHalf = sqrt(abs(D))

sgnp = -1 when ((-b/2.0) + dHalf) < 0, 1 otherwise    
sgnq = -1 when ((-b/2.0) - dHalf) < 0, 1 otherwise

p = 0.0 when ((-b/2.0) + dHalf) = 0,
     sgnp*(abs((-b/2.0) + dHalf))^(1.0/3.0) otherwise
q = 0.0 when ((-b/2.0) - dHalf) = 0,
     sgnq*(abs((-b/2.0) - dHalf))^(1.0/3.0) otherwise

s = (-b/2.0)/sqrt(-a*a*a/27.0);  theta = acos(s)

-- roots of the form x + iy

x1 = 2.0*p - (c1/3.0) when numRealRoots = 1 and abs(D) < 1.0e-10,
     (p+q) - (c1/3.0) when numRealRoots = 1 and abs(D) > 1.0e-10,
     2.0*sqrt(-a/3.0)*cos(theta/3.0) - (c1/3.0) otherwise

x2 = -p - (c1/3.0) when numRealRoots = 1 and abs(D) < 1.0e-10,
     -(p+q)/2.0 - (c1/3.0) when numRealRoots = 1 and abs(D) > 1.0e-10,
     2.0*sqrt(-a/3.0)*cos((theta/3.0) + 120) - c1/3.0 otherwise 

x3 = x2 when numRealRoots = 1,      
     2.0*sqrt(-a/3.0)*cos((theta/3.0) + 240) - c1/3.0 otherwise 

y1 = 0.0  -- no matter what, must have at least one real root

y2 = (p-q)*sqrt(3.0)/2.0 when numRealRoots = 1 and abs(D) > 1.0e-10,
     0.0 otherwise

y3 = -y2 when numRealRoots = 1 and abs(D) > 1.0e-10,
     0.0 otherwise

root1 := {x1,y1}:;     root2 := {x2,y2}:;   root3 := {x3,y3}:

------- Enter the coefficients for the cubic here --------------
--            a0*x^3 + a1*x^2 + a2*x + a3 = 0

a0 = 35;    a1 = 15;    a2 = -5;    a3 = -20

root1:{0.76,0.00}
root2:{-0.59,0.64}
root3:{-0.59,-0.64}

----------- Check the solution------------
include "::incl:complex ops"
z(x) = a0*Ccube(x) + a1*Csqr(x) + a2*x + {a3,0} -- The complex cubic

-- confirm that z(x) is zero at the roots
z(root1):{0.00,0.00}
z(root2):{0.00,0.00}
z(root3):{0.00,0.00}

-------- check the solution graphically
-- define an array that samples points on the surface of:
--    Z = abs(z(x)) vs. X = real(x) , Y = imaginary(x)
-- This surface should dip to zero at the roots. (The sampled surface will dip near zero but in general is not sampled exactly at a root).

surf[ix,iy] =  x:=Cscale(ix,iy), Cabs(z(x)) dim[m,m]

-- The index for surf[] runs from 1 to m. Scale the index to get real and imaginary parts ranging from Xmin to Xmax and Ymin to Ymax
scale(i,lo,hi,nsteps) = (i-.5)*(hi-lo)/nsteps+lo
Cscale(ix,iy) = {scale(ix,Xmin,Xmax,m), scale(iy,Ymin,Ymax,m)}

m=80;   -- surface is sampled on an m by m grid
Xmin = -1; Xmax = 1
Ymin = -1; Ymax = 1
Zmin =  0; Zmax = 50
image surf                  -- show image of the surface
plot {root1,root2,root3}    -- show locations of roots

plot z({X,0})[1]/Zmax     -- plot the real part of z for real x
-- This should be zero at real roots. On the plotted surface, real roots are located along y=0 so the real cubic plotted in this way should pass though its real roots.
plot 0


